ON THE EFFECT OF GHOST FORCE IN THE QUASICONTINUUM 
METHOD: DYNAMIC PROBLEMS IN ONE DIMENSION 
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Abstract. Numerical error induced by the "ghost forces" in the quasicontinuum method is 
studied in the context of dynamic problems. The error in the W 1 ' 00 norm is analyzed for the time 
scale O(e) and the time scale O(l) with e being the lattice spacing. 

Key words. Quasicontinuum methods; Ghost force; Dynamic problems. 

AMS subject classifications. 65N15; 74G15; 70E55. 

1. Introduction. The present paper is mainly concerned with the error pro- 
duced by the ghost forces in quasicontinuum (QC) type of multiscale coupling meth- 
ods for crystalline solids. In these methods, one reduces the degrees of freedom of an 
atomic level description by replacing part of the system with continuum mechanics 
models (31] |2j [T5J EHl ED] ■ Such integrated methods have been very useful in studying 
mechanical properties of lattice defects. It allows one to simulate a relatively large 
system while still able to keep the atomistic description around critical areas, such as 
crack tips and dislocations cores. These methods have also drawn a great deal of at- 
tention from numerical analysts. We refer to JTSJ [TTJ [TUJ [35J 03 [HI [Ml [23] and references 
therein for a partial list of the representative work. Nevertheless, many challenges in 
the analysis of these methods still remain. Examples include high-dimensional prob- 
lems, systems with line or wall defects, and solutions near bifurcation points. We 
refer to [HI [26) for a review of the state-of-art of this field. A critical issue that arises 
in the numerical analysis is the ghost force, which is the non-zero forces on the atoms 
near the interface at the equilibrium state |30j . For statics problems, the elimination 
of ghost forces has been a necessary ingredient to achieve uniform accuracy [TDJ [25] . 

In the static case, the QC method couples a molecular statics model to the 
Cauchy-Born elasticity model. For one-dimensional models, the influence of ghost 
forces has been explicitly characterized in [2TJ [551 E| • They found that the ghost force 
may lead to an 0(1) error for the gradient of the solution, and the width of the result- 
ing interfacial layer is of the size C(e|lne|), where e is the equilibrium bond length. 
The influence of the ghost force for a two-dimensional model with planar interface 
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has recently been studied in 0]. It was found that the ghost forces still lead to an 
0(1) error for the gradient of the solution, while the interfacial layer caused by the 
ghost force is of the size 0(y/e), which is much wider than that of the one-dimensional 
problems. 

The QC method can be extended to dynamic problems using the coarse-grained 
energy and the Hamilton's principle [5^1 H7]. The dynamic QC method couples an 
elastodynamics model with a molecular dynamics model. Many dynamic coupling 
methods with similar goals have been developed JJ [SJ [5J [7J [HJ [T71 [351 [33J 133 [35] ever 
since. However, very little has been done to address the stability and accuracy of these 
methods. Most numerical studies have been focused on the artificial reflections at the 
interface. The reflection is caused by the drastic change in the dispersion relation 
across the interface, which is often due to the difference between the mesh size in 
the continuum region and the lattice spacing in the molecular dynamics model. The 
reflection can be studied by considering an incident wave packet traveling toward the 
interface and examine the amplitude of the reflected waves. The issue of ghost forces, 
however, has not yet been addressed. 

The purpose of this paper is to study the effect of ghost forces in the context of 
dynamic problems. Motivated by the results in the static case, we expect that ghost 
forces will continue to play an important role in dynamic coupling models. To focus 
primarily on the issue of ghost forces, we consider the dynamic model j29l[27] derived 
from the original QC method when the mesh size coincides with the lattice spacing. 
In addition, the initial displacement is given by a uniform deformation. This allows 
us to compute the error caused only by the ghost forces. The error will be studied 
in the W 1 '°°-norm as in the static problem [211 E E21 H]- Our study shows that the 
error, which is initially zero, grows very quickly, and already becomes 0(1) at the 
time scale 0(e) . The error exhibits fast oscillations, with amplitude on the order of 
e. On the time scale 0(1), which is typically the time scale of interest, the amplitude 
of the oscillations grows, and it is bounded by an 0(y/e) quantity. The average of the 
oscillations has a peak at the interface. In contrast to the static case, where the error 
is mainly concentrated at the interface, the error in the dynamic case is observed in 
the entire domain on the time scale 0(1). These observations are quite different from 
those of wave reflections, and it indicates that the effect of ghost forces is a separate 
numerical issue. 

The rest of the paper is organized as follows. In Section 2 we describe the one- 
dimensional atomistic model and the derivation of the QC model, and briefly demon- 
strate the appearance of the ghost forces. In Section 3 we show results from several 
numerical tests. They provide some insight into the evolution of the error. The next 
three sections are devoted to the analysis of the error for short and long time scales. 
We draw some conclusions in the last section. 

2. Motivation and the Formulation of the Problem. As in [TJ] , we consider 
the dynamic problem of a one-dimensional chain of atoms. The interatomic interaction 
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is assumed to be among the nearest and the next nearest neighbors. Let x be the 
reference position of an atom, and y(x,t) be the current position at time t. The 
equations of motion for the atoms in the chain read, 



y(x,t) ~ C at [y\(x,t) = 0, ieL 
y(x, 0) = x, y(x, 0) = x, 
y(x, t) — x is periodic with period 1. 



(2.1) 



Here, we have set the mass to unity, and L = {je,j eN}fl (—1/2, 1/2) with e being 
the lattice parameter. The operator £ at is defined as 



£ at [z](x,i) 



1 [k2z(x — 2e, t) + Kiz(x — e,t) + k.\z{x + £,€) + K2z(x + 2e, t) 



2(ki + K 2 )z(x,t)]. 



(2.2) 



Since the issue of ghost force arises even for harmonic interaction, we consider 
here a linear model, which can be considered as a harmonic approximation of a fully 
nonlinear model. In (|2.2[) . k\ and k 2 are the force constants computed from an 
interatomic potential. For example, for a pair potential, the energy is given by 



E 



E 



9 



y(x + s) -y{x) 



y(x + 2s) - y(x) 



Direct calculation yields 

Ki = ^"(1), n 2 = <p"(2). 

One commonly used model is the Lennard- Jones potential [16j : 

ip(r) = (cr/r) 12 - {a/rf . 

If only the nearest and the next nearest neighborhood interactions are considered, the 
lattice parameter is given by 

-12 \ 1/6 



2 l/6 



1 



1 



In this case, the force constants are 

Ki = 156C 2 - 42C w 18.886 and n 2 = 2~ 6 (156C 2 2~ 6 - 42C) 



-0.323, 



where C = (l + 2 _6 )/[2(l + 2 -12 )]. The above formula has also appeared in [13] . 
Notice that the second force constant is negative, but it is much smaller than the first 
force constant in magnitude. 

With the harmonic approximation, the potential energy takes the form of 

2 / / , n \ r \ \ 2" 



M ( y{x + e) - y(x) 

2 



«2 (y{x + 2e) - y(x) 
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The dynamic model (|2.1[) can be derived from this energy using Hamilton's principle. 
Notice that the energy can be divided into the energy at each atom site: i.e. E = 
^2 x E(x), in which 



K2 fy{x + 2e) - y(x)\ 2 m f y(x + e) - y(x) x 2 



K i ( V( x ) ~y(x~£)\ 2 , k 2 ( y(x) - y(x - 2e) x 2 



E(x) 



In the QC method, one defines a local region where the atomistic model is approx- 
imated by the Cauchy-Born elasticity model [3]. One also defines a nonlocal region 
where the atomistic description is kept. Without loss of generality, we assume that 
the interface is located at x = and the nonlocal region is in the domain x < 0. We 
further assume that the mesh size is equal to the lattice parameter to primarily focus 
on the effect of ghost forces. The Cauchy-Born approximation of the energy in the 
local region is given by 

ki+4k 2 f y(x + e) - y(x)\ 2 m + 4k 2 f y(x) - y(x - e) x 
= — 4— ( S J +^ 

At the interface x — 0, the energy takes a mixed form: 
E(0) = K -l 



K2 ( y(x) - y(x - 2e) \ 2 Ki ( y(x) - y(x - e) x 2 



rei + 4k 2 / + e) - j/(aj) x 2 



With such energy summation rule, we may write the QC approximation of £ at as 
£ qc , with £ qc given below. For x < — 2e, 

£ qc [z](a;, t) = f -2 [k22(.t — 2e, t) + n\z(x — e,t) + k\z{x + e, t) + K2z(x + 2e, t) 

- 2(«i + K 2 )z(x,t)] , 

and for x > 2e, 

C qc [z](x,t) = C C b, Ccb = £ _2 (ki +4:K 2 )[z(x - e,t) - 2z(x,t) + z(x + e,t)}. 

This is exactly the operator corresponding to the Cauchy-Born approximation. 
At the interface, we have for x = — e, 

C qc [z](x, t) = e~ 2 [k 2 z(x — 2e, t) + K\z(x — e,t) + K\z{x + e,t) + ~^ z i x + 2e, t) 

- (2K 1 +3K 2 /2)z{x,t)], 



for x = 0, 



£qc[z](x,t) — e 2 K2z{x — 2e, t) + K\z{x — e, t) + (ki + 4k2)z(x + e, t) 
- (2ki + 5k 2 )z(x 7 t) 
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and for x = e, 
C qc [z] (x,t) 



= e- 2 



— z[x — 2s, t) + (fti + Ak 2 )z(x — e, t) + (ki + Ak 2 )z(x + e, t) 



- (2/tx + I7n 2 /2)z{x,t) . 
Using the Hamilton's principle, we write the QC model as 

y(x,t) - C qc [y](x,t) = 0, x e L, 
0) = x, y(x,0) = x, 

y(x, t) — x is periodic with period 1. 



(2.3) 



The initial and boundary conditions have been chosen as a uniform deformation 
in order to identify the effect of the ghost force. We will compute the deviation 
of the solution away from the equilibrium. For this purpose, we define the error 
y(x, t) = y(x, t) — x, and we have, 



y{x,t) - C qc [y](x,t) = y(x,t) - C qc [y - x](x,t) 
= C qc [x](x,t) = f(x,t), 



(2.4) 



with / given explicitly by 



if |x| > 2e, 



f(x,t) = { 



_K2 

e 

2k 2 

e 

l<2 



if x = — 
if x = 
if x = e 



e, 



(2.5) 



The function /(x, t) is precisely the ghost force. Since it is independent of the temporal 
variable, we denote it by f(x) for simplicity. Finally, we supplement the above problem 
with the homogeneous initial condition and periodic boundary condition as 



y(x, 0) = and y(x, 0) = 0, x 6 L. 
y(x,t) is periodic with period 1. 



(2.6) 



3. Observations from Numerical Results. Since the operator £ qc coincides 
with Cqb in the local region, and with £ at in the nonlocal region, it is natural to look 
at models similar to (|2.4[) , in which C qc is replaced by either Cqb, or C a t in the entire 
domain. Therefore, our numerical experiments are conducted for the following three 
models: 

• Model I. C qc is approximated by £cb: jj — £cB[y] = /. 

• Model II. C qc is approximated by £ at : y — £ a t[y] = /■ 

• Model III: The quasicontinuum model (|2.4p . 

In all these models, we impose homogeneous initial condition and periodic boundary 
condition (J2.0D . 
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As an example, the force constants are obtained from the Morse potential [23]. In 
particular, we choose K\ = 4.4753 and k 2 = 0.4142. All the simulations are performed 
in the domain x e [—1/2, 1 /2] , and the ODEs are integrated using the Verlet's method. 
Since all three models are Hamiltonian systems, this method is particularly suitable. 

We first show the solutions computed from the three models at different time 
step. The results are shown in Fig. 13.11 For this set of numerical tests, we have 
chosen e = 1/2000. We observe that the error first developed at the interface, and 
then it starts to spread toward the local and nonlocal region for all three models. 
Another noticeable feature is that the error exhibits a peak at the interface, and the 
peak remains for all later time. At t = 1, the error is observed in the entire domain. 
Our main observations here can be summarized as following: (1) In the presence of 
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Fig. 3.1. The gradient of the error. Left to right: Solution computed from Model I, II and III. 
From top to bottom: The solutions at time t = 0.01,0.05, 0.2 and 1. 



ghost forces, the error grows very quickly. It reaches 0(1) on the time scale of O(e); 
(2) At the time scales of t = 0{e) and t = 0(1), the solutions of all three models are 
qualitatively the same. 

Next we monitor the solution for those atoms near the interface. In Fig. 13. 2\ we 
show the time history for those atoms. We observe that for most of the time, the 
error oscillates around certain constant values, and the constant values depend on the 
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location of the atom. These constant values show a peak at the interface x = 0. 



i T 
i r 






Fig. 3.2. T/ie time history of the gradient of the error near the interface. Left to right: 
Solution computed from Model I, II and III. From top to bottom: The solutions near the interface: 
x = — 3e,x = —2e,x = —e,x = 0,x = e and x = 2e. 



In the last two figures, Fig. 13.31 and Fig. 13.41 we show the time history of the 
solution at the interface for various values of e. The main observation is that the 
amplitude of the oscillation decrease as e gets small. However, the constant values 
around which the solutions oscillate do not change as e varies. 

4. Explicit Solutions for the Approximating Model. In view of the nu- 
merical results, it seems that the solution of model I bears similarity to the dynamical 
behavior of the original problem (|2.4[) on the time scale 0(e) and time scale 0(1). 
Therefore we will turn to this model to study the effect of ghost forces. Model I is 
convenient to analyze, particularly because it admits explicit solutions of a simple 
form. Recall that in Model I, we solve the following problem, 

y(x,t) - £cb[v] (x, t) = f(x), x 6 L, 

y(x,0) =0,y(x,0) = 0, (4.1) 

y(x,t) is periodic with period 1. 
7 




Fig. 3.3. The gradient of the error. The solutions are computed from Model I with different 
choice of e. Left to right: Solution computed for e = 1/2000, e = 1/4000 and e = 1/8000. From top 
to bottom: The solutions near the interface: x = — 3e, x = — 2e, x = — e, x = 0, x = e and x = 2e. 



Without loss of generality, we let L = (—1/2, 1/2] with N atoms. We assume that N 
is an even integer for technical simplicity. Obviously, e = 1/N. We will switch to the 
notation that, 

w(n,t) =w(-l/2 + ne,t), n = 1, 2, ■ ■ ■ , N. (4.2) 

We now express the solution of (|4.1[) in an explicit form. To begin with, we 
consider the lattice Green's function, which is defined as the solution of the following 
problem: 

' G{n,t)-C C B[G\(n,t) - 0, n = 0,...,N, 

G(n,0) = 0,G(n,0) = 5 n , (4.3) 
G(n, t) = G(n + N, t). 

Given this Green's function, the solution of (|4.1j) is given by 
y(n, t)= V G ( n ~ m ' 1 ~ s )/( m ) ds ' 

J Vm=l / 
8 





Fig. 3.4. The gradient of the error. The solutions are computed from Model III with different 
choices of e. Left to right: Solution computed for e = 1/2000, e = 1/4000 and e = 1/8000. From 
top to bottom: The solutions near the interface: x = —3e, x = —Is, x = —e, x = 0, x = e and 
x = 2e. 



where / is given by (|2.5[) under the transform (|4.2I) . As a result, we have 

JV 

V G{n - m, t)f(m) = — (2G(n - L,t) - G(n- L - 1, t) - Gin - L + 1, t)) , 
1 e 

where we have set L = N/2. 

By separation of variables, we have the following explicit form for the lattice 
Green's function G: 

. t 1 v-^ 1 sinfwttl 2knir ,, 
Gin, t) = 1 > — cos , (4.4) 

fc=i 

with LUk being the dispersion relation given by uik = (2/e)y / /ti + 4^2 sin(fc7r/iV). 
Using X) m =x ^( TO ) = °' we write 

E , . , . 4k 2 v 1 ^ 1 sin[w fc (i- a)] . 2 kir 2kn 
G{n-m,t- s)f(m) = — ^ — sin — cos — (n-L). 

m=l fe=l 

9 



This leads to 

/ \ /"* f s-n , ,, ,\ , e 2«2 v-^ 1 . o Wfet 2fc7r 
y(n,t) = / > G(n — m, t — s)j[m) us = > sin cos (n—L). 

Using the above expression we bound y(n,t) as 

\y(n,t)\<-^-e. (4.5) 

Ki + 4« 2 

This estimate shows that the magnitude of the error y(n, t) is as small as 0(e) for all n 
and all time t. This in turn suggests that the error induced by the ghost force is small 
in the maximum norm, which is consistent with that of the static problem [5J 1221 [3] . 
Next we consider the discrete gradient of the error. A direct calculation gives 

n , .s y(n + l,t) - y{n,t) 
Dy(n,t) = 

£ 

N—l 

1 4/-C2 ^ . 2kir kir . 2 u k t 

= — — - — > sin — — (n + 1/2 — L) sin — sin — — . 

N ki + 4k 2 ^ N v 7 ' AT 2 
fc— l 

Clearly we may write the above expression as 

, , 1 Ak 2 x— > . 2fc7r . . TN . kir . , Wfci ., „. 

Dy(n, i) = — > sin (n + 1/2 - L) sin — s n 2 — . 4.6 

1 z k=0 

It follows from the above expression that Dy(n,t) is anti-symmetric in the sense of 

Dy(n, t) = —Dy(N - n - 1). 
Therefore, we only consider the case n > L. By (|4.5p we bound Dy(n,t) trivially: 

4|«a| 



|£y(M)l < 



rti + 4k 2 



This shows that Dy(n,t) is uniformly bounded for all n and all time t. In the next 
two sections, we seek for a refined pointwise estimate of Dy(n, t) when t is of 0(1) and 
of 0(e). Notice that the same method can be employed to obtain a refined pointwise 
estimate of y(n,t). We leave it to the interested readers. 

5. Estimate of the error over long time. In this section, we estimate the 
error for t = 0(1). By (|4.6p , we write Dy(L,t) as 

, . 1 4kt x— -v 9 kir 9 uji.t 

Dy(L,t) = - — -I—} sin 2 — sin 2 

^ ' ; TVki +4k 2 4^ ^ 2 



Using the identity 



2 kir _ N 



N 2 

fe=0 
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Wc write 



Dy(LA) = — — - — > sin — — + — — > sin — cos(Wfci) 

yv ' ' TV m + 4k 2 ^ TV TV ki + 4k 2 ^ TV v ' 

k—0 A;— 



ki + 4k 2 TV ki + 4k 2 ^ TV 

k— 



When n ^ L, we use the fact that 

^ sin ( n + 1/2 _ L) sin — = 0, 



TV v ' y TV 
and we write the expression of Dy(n, t) in (|4.6p as 

1 2^2 2fc7T /C7T 

^MM) = T7 — — j— Z^ sm ^r( n + 1 / 2 ~ £)sm — cos(Wfci), 
TV /ci + 4k2 f-^ TV TV 

fc— o 

which can be further decomposed into 

Dy(n, t) = -—^-Y, ^Jf{^ + — (» + 1/2 - L) 

k=0 

- sin L k t - ^(n+ 1/2 -lM. (5.2) 
To bound Dy(n, t), we need to estimate an exponential sum of the following form, 

N 
k=0 

where the shorthand e(/(fc)) = exp(27ri/(fc)) is assumed. The basic tool that will 
be used is the truncated form of the Poisson summation formula due to Van der 
Corput |34j . The following form with an explicit estimate for the remainder term can 
be found in P3J Lemma 7]. 

Theorem 5.1. (Truncated Poisson) Let f and <fi be real-valued functions 
satisfying the following conditions on a closed interval [a, b] : 

1. f" and 4>'(x) are continuous; 

2. < f[x) < C ; 

3. There are positive constants H, U,4>q,4>\, \ such that U > 1, < b — a < ACT 
and 

\4>ix)\ < 4> H, \4>'ix)\ < faH/U. 
For any A,0 < A < 1, the equation 

£ #n)e(/(n)) = £ / </>(x)e(f(x)-ux)dx + 0R (5.3) 

a<n<b a-A<u<P+A^ a 
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holds, where a = f (a) , (3 = /'(&) and 

R= (0 O + A0 1 )i7^9.42 + 9C O + 12A + ^(lOA^ 1 + 2 In A" 1 

+ 4.5(1 + A)" 1 - 4.5 ln(l + A) + 6.5 ln(/3 - a + 2)) 

Here 9 is a function such that \9\ < 1. 

The assumption /" > can be relaxed to either /" > or /" < 0. In the latter 
case, the second condition is replaced by —Co < f"{x) < 0. 

5.1. The estimate for Dy(N/2,t). To bound Dy(N/2,t), we start with ([53]) . 
Based on the above theorem, we transform the exponential sum in (|5.ip to a shorter 
sum with a bounded remainder. To clarify the depcndance of the constant, we denote 



7 = t\/ni + 4k 2 , 

and assume that 1 < 7 < N since t = O(l) and N > 2. We also denote by [p \ the 
integer part of a real number p, and denote its fractional part by {p} = p — \_p \. 
Lemma 5.2. Let cf>(x) = sm 2 (irx/N) and f(x) = j/{ire) sm(irx/N). There holds 



1 

N 



N 



E<^m/(*o) 



k=0 



< 



2 



E 



q-1/2<i/<7+1/2 



N/2 



(x)e(f(x) — vx) dx 



+ Ce(l + £ 7 + log( 7 + 2)), 

where C is independent of N,t and 7. 

Proof. It is easy to write the exponential sum into 



(5.4) 



N 



N/2 



E <K fc M/(fc)) = 2 E m<f(k)) - 4>{N/2)e{f{N/2)). 

fe=0 fc=0 

With such choice of / and 0, we have 

a = 0,b = N/2,a = 0,f3 = 'y, 
00 = 01 = 1, # = 7T, A = 1, f7 = JV/2. 

Setting A = 1/2 and using Theorem 15. 11 we obtain 

W/2 -, ,N/2 



fc=0 q-1/2<i/<7+1/2 -/U 



(x)e(f(x) -ux)dx + 9— , 



where 



i? = 2tt(9.42 + 9?r7/iV + 6 + tt" 1 (20 + 2 In 2 + 3 - 4.5 ln(3/2) + 6.5 ln( 7 + 2) 
This immediately implies that there exists a constant C such that 

|<9i?/iV| < Ce (1 + e 7 + 111(7 + 2)) • 
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We obtain (|5.4[) by combining the above two inequalities. □ 

Remark 5.3. The choice of A is not unique. However, it cannot be too small. 
Otherwise, the remainder term blows up as A — > 0. 

To bound the shorter sum in (|5.4[) . we shall rely on the following first derivative 

test. 

Lemma 5.4 (First derivative test). |Hl Lemma 1, p. ^7/ Let r(x) and d(x) 
be real-valued functions on [a,b] such that r(x) and 0'(x) are continuous. Suppose 
that 9'{x)/r{x) is positive and monotonically increasing in this interval. If < Ai < 
9'(a)/r(a), then 



r(x)e(6(x)) Ax 



< 



1 

ir\i 



Remark 5.5. If9'{x)/r(x) is negative and monotonically decreasing on [a,b] and 

e'(a)/r(a) < Ai < 0, 

then we obtain the same bound by taking complex conjugates. Moreover, if 8' (x)/r(x) 
is monotone on [a, b] and 

\6'{x)/r{x)\ >X X >0 x€(a,b), 

we obtain the same bound by combining the above two cases. 
We write 



1 



N/2 



4>{x)e{f(x) — vx) dx 



1 



N Jo TT Jo 



r/2 



(p{y)e{G v {y)) dy, 



where (p(y) — sin 2 y, and G„(y) — ( N/n) (7 sin y — vy) for any v £ N. We define 
Fy{y) = G' u (y)/ip(y). By Lcmma [5~2l it remains to estimate the integral J^ 2 (p(y)e(G v (y)) dy 
for v = 0,-- • , |_7 + 1/2J. The three cases v = 0, v = 1, • • • , [7 + 1/2 J - 1 and 
v = L7 + I/2J will be treated separately in the following lemmas. 
Lemma 5.6. There holds 



tt/2 



<p{y)e (G {y)) dy 



< 2{N 1 )- 1 ' 2 



(5.5) 



Proof. For any 8 € (0,7r/2) that will be determined later on, we have, for any 
y€(0,n/2-5), 

F (y) > F {ir 2-6) = —L — ^ -^tan<5> — !-, 

7T COS z IT 7T 

where we have used the fact that tanx > x for x € [0,7r/2]. Using Lemma 15.41 with 
Ai = N^8/-k, we obtain 



iv/2-8 



<p{y)e (G (y)) dy 



< 



1 
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The integral over the complementary portion of the interval can be bounded trivially: 



ir/2 



ip{y)e (G (j/)) dy 



Itt/2-8 

On adding the two estimates we deduce that 

r/2 

<p(y)e (G (y)) dy 



• it/2 

I dy<5. 

tt/2-S 



< 



Nj5 



This is minimized by taking 5 = {N^)^ 1 / 2 <E (0,7r/2), and we obtain (|5.5I) . □ 
The second case is more involved since F v changes sign over (0,7r/2). 
Lemma 5.7. Ifl < v < 7, tften 



ix 12 



<p{y)e{Gv{y))dy 



< 



3tt 

//V7' 



(5.6) 



Proof. For 1 < v < 7, there exists y v € (0, tt/2) such that F v (y v ) = with 
cosy„ = 1//7. For any 77 S (0, min(j/,/, 7r/2 — y v )) that will be chosen later, we write 



7r/2 



<p(y)e (G v (y)) dy 



v{y)e{G v {y)) dy. 



We deal with the three integrals separately. 

Using Lemma 15 .41 with Ai = \F v (y v — rj)\, we obtain 



v{y)e(G v {y)) dy 



< 



and 



\F v {y v -r,)\ = F v {y v -r,) = 



Tr\F„(y u - 77)|' 
Nj cos(y v - 77) - cosy,, 



> 



> 

tt z sin y v ' 

where we have used Jordan's inequality 



sinx > —x 

7T 



7r sin 2 {y v - 77) 
2iV 7 sin(y„ - 77/2) sin(?7/2) 

7r sin 2 (y v - rj) 

2./V7 sin(77/2) 
7r sin y v 



x e [0,tt/2]. 



This gives 



Vv—v 



<p{y)e(Gv(y)) dy 
14 



< 



7r sin y v 
~2Nyq' 



(5.7) 



Using Lemma again with Ai = \F v (y v +T))\, we have, for the second integral, 



it/2 



ip(y)e(G v (y)) dy 



< 



1 



Furthermore, 



\F»(Vv + V)\ = -FAVv+v) 



*\Mv» + v)\ 



N'y cos y v - cos(y„ + 77) 



> 



7r sin 2 (2^ + rf) 
2Nj sin(y v + T]/2) sm(r]/2) 

7T sm 2 (y v +r)) 
2Njr] sin(y v + 77/2) 
7r 2 sm 2 (y^ + 7^) 



This leads to 



it/2 



<f(y)e(G v (y)) dy 



< 



7T sin 2 (7^ + 77) 
2^777 sin^ + 77/2) ' 



If y v £ (0,7r/4], we would require that 77 e (0,y v )- We will have, sm(y v + 77) < 
aia2y v < 2smy L/ and siny^ < siii(?/„ + 77/2) since y u < y v + rj/2 < 2y v < w/2. The 
bound for the above integral is changed to 



it/2 



ip(y)e(G u (y)) dy 
We estimate the remaining integral trivially: 



< 



ip(y)e{G v (y)) dy 



2ir sin y v 
N^/rj 



< 277. 



Vu-V 



Summing up all the above estimates, we obtain 



it/ 2 



*p(y)e(G v (yj) dy 



5tt sin y u 

< — + 277. 

- 2^777 ' 



Taking 77 = (N'y) siny u , which is less than y v , we obtain 



r/2 



ip(y)e(G v (y)) dy 



< 



• r )7T 



2V^7 y/Nj 



2 sin yv ^ 5tt/2 + ^/2 3tt 



On the other hand, if y v E (7r/4,7r/2], wc would require that 77 <E (0, 7r/2 — y v ). 
We will have sm.y v < sm(y v + 77/2) since y v < y v + rj/2 < y v + 77 < 7r/2, and 
sin 2 (y I y + 77) < 1 < 2 sin 2 y v since sin 2 y v > 1/2. We bound the second integral as 



tt/2 



<p(y)e (G„(y)) dy 



Vv+V 



< 



7r sin y v 
Njrj 
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This yields 



tt/2 



<p(y)e{G v (y)) dy 



3nsiny u 

< h 2n. 

~ 2Njtj ' 



In this case, we can choose rj = (Nj) 1//2 provided that 

77 < tt/2 — y v . 

This immediately implies that 

r/2 



(5.8) 



¥>(y)e(G„(y)) dy 



< 



Such choice of 7? is feasible since 



v > 1 > v/tM 

which yields 1/ > 777, or equivalcntly, 77 < cosy„ = sin(7r/2 — y„) < tt/2 — y v . This 
directly gives (|5.8|) . Finally we get (|5.6[) . □ 

Next we consider the endpoint case v = [7 + 1/2 J. 

Lemma 5.8. Let v = [7 + 1/2 J . Ifv>j, then 



tt/2 



<p(y)e(G y (y))dy 



< 



If v < 7, then 



r/2 



ip{y)e{G v {y)) dy 



< 



3tt 

//V7' 



(5.9) 



(5.10) 



Proof. If ^ > 7, then there exists a stationary point y„ of F v (y) with cosy„ = 
v/j — \/v 2 j^f 1 — 1. Because F v {y) is monotonically increasing over (0, y„) and mono- 
tonically decreasing over (y„,7r/2), we get min yg ( 0jW /2) |^(y)| > |-FMz/y)l- To each of 
the two intervals we apply Lemma IST41 with A = \F v (y v )\. On adding these estimates 
we deduce that 



tt/2 



which yields (|5.9I) by using 



ip(y)e(G v (y))dy 



< 



(5.11) 



= ^^/7+^ 2 /7 2 -l)>^. 



When 1/ = L7 + 1/2J < 7, we invoke the estimate ff5~6]) to get (|5.10jl . □ 
Summing up the above estimates for the shorter sum, we obtain the estimate for 
the exponcntinal sum in (|5.4p . 
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Lemma 5.9. There holds 

N 



1 

N 



5>(fc)e(/(fc)) 



<C^(l + £ 7 + log( 7 + 2)) + ^ + ^ , (5.12) 



where C is independent of N,t and 7. 

Proof. Denote vq = [7 + 1/2 J. If Vq > 7, then 



"a pn/2 fir/2 "0-1 

2 / <Piv)e {Gu{y)) dy = / <p(y)e (G (y)) dy + £ / 
^o- 70 ■'o „ =1 Jo 



(p(y)e(G u (y)) dy 



tt/2 



v(y)e{G Vo {y)) dy. 



Using the estimates (|5.5[) , (|5.6[) , and (|5 . 9[) , we bound the right hand side of the above 
sum as follows, 



"a „ir/2 

V/ <p{y)e{G v {y)) dy 
„- n Jo 



i/Q-l 



< 2 (iV 7 )" 1/2 + 3tt ( N l)~ 1/2 + 4(^7)-! 

< 7r(7V7)- 1/2 (3 7 + 1) +A(N-f)- 1 . (5.13) 



If uq < 7 , we have fo = L 7 J ; an d 



£ / p(y)e (G„(y)) dy = / <p(y)e (G (y)) dy + £ / 
Proceeding along the same line that leads to (|5.13[) . we obtain 



it/2 



V(y)e (G v (y)) dy. 



J2 <p(y)e(G v (y)) dy 



< 2(A^ 7 )" 1/2 + 37r( 7 /A0 1/2 . (5.14) 



Combining the estimates (|5.13|) . (|5.14|) and (|5.4j) . we obtain f)5.12[) . □ 

Substituting the estimates (|5.12[) and (|5.4p into (|5.1j) . we obtain the pointwisc 

estimate for D(N/2, t) as follows. 

THEOREM 5.10. J/i = 0(1) and n = N/2 or N/2 - 1, then 



Dy(n,t) 



H 2 



Kl + 4«2 



< C- 



«2 



Kl + 4«2 



e [l + £ 7 + log(7 + 2)] + A /-(l + 7) ), (5.15) 



where C is independent of N,t and 7. 

The above estimate means that Dy(n,t) is in the O(-Je)— neighborhood of 
K2/(«i + 4k 2 ) when n = N/2 orn = A^/2 - 1. 

5.2. The estimate for Dy(n,t) with n ^ L. By (|5.2j) . we need to estimate 
two exponential sums X^fclo ^(^OK/CO) with 

. 7tx 7 , n ji + l/2-L 
ffl x = sm — , tlx) — — sm x, 
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and 



. TTX ., . 7 . TTX 11+1/2- L 

mix = sm — , fa; = — sm x. 

VK ' N' JK ' ire N N 

In what follows we only give the full details for estimating of the first exponential 
sum, and the same proof works for the second exponential sum. Denote by g = 
(n + 1/2 — L)/N, proceeding along the same line that leads to (15.41) and choosing 
A = max(l/2, 1 — { 7 + g}), we get 



1 

N 



N 



k=0 



< 



+ Ce(l + £7 + log( 7 + 2)), 



(5.16) 



where C is independent of N,t and 7. Here tp(y) = smy and G v (y) = ( N/ir) (7 sin y - 
gy - vy). We also define F v {y) = G' v (y) / ip{y) . 
Lemma 5.11. There holds 



zr/2 



<p(y)e(G (y)) dy 



<nun(2 ( AT7)-^- TT L- x ). (5, 



17) 



Proof. If 2(^7) x / 2 < l/(n + 1/2 — L), then we proceed along the same line that 
leads to Lemma l5~6l to obtain 



it/2 



v(y)e(Go{y)) dy 



< 2{N 1 )- 1 / 2 . 



Otherwise, using Lemma 15.41 with A = (n + 1/2 — L)/tt, wc obtain 



tt/2 



<f(y)e(G {y)) dy 



< 



n + 1/2 - L 



Combining the above two inequalities gives (|5.17|) . □ 
Lemma 5.12. Let v = [7 + g\ +1, then 



tt/2 



tp(y)e{G v (y)) dy 



< A(Nj)~ 2/3 . 



(5.18) 



Proof. The function F v has a stationary point y v that satisfies cosy^ = ^ j(y — g). 
In this case, applying the first derivative test directly to the integral may yield a 
bound of the form 1/(N siny„), which is undesirable since y v can be very close to 
zero when v is close to 7 + g. Instead, for any 8 G (0, tt/2) to be determined later on, 
we have 

(\{y)e(G v {y))dy 
Jo 

18 



< 



<p{y) dy < 



j s 2 
ydy = —■ 



If Vv < S, then we use Lemma with A = |F„((5)|. This gives 



7T/2 



v{y)e{G v {y))dy 



< 



1 



If j/„ > <5, then we proceed along the same line that leads to (|5.11[) to get 



Tl/2 



< 



TTmm ye [s t7r / 2 ]\F v (y)\ w\F v (y„)\ 



A direct calculation gives 
N 



\MS)\ 



7T 

N 



, 1 — cos y v cos S N 1 — cos <5 

- e) :— 7 > — - q) — r-j— 

sin o ir sin o 



(f — g) tan — > 



5 ^_ N{v- e ) 



2tt 



5, 



and for y u > 5, 



N 



\F v {y v )\ = —{y- o)smy v > 



2N(v - g) 2N{v - g) 
=5 V" > Z5 ' 



Combining the above four inequalities, we obtain, for any S <G (0, 7r/2), 



tt/2 



ip(y)e(G v {y))dy 



< max 



ttI^^I'ttI^^)!; - N(v-q)6 



< 



To sum up, we have 



tt/2 



tp(y)e(G u {y)) dy 



S 2 
< — 



2 N{v-q)5 
Taking 5 = ^^(N^ - g)) _1/3 € (0,tt/2), we get 



tt/2 



tp(y)e(G„(y))dy 



<^ /3 [^-fi)r !/3 <4[%-r J/3 > 



(5.19) 



(5.20) 



which yields (|5 . 1 8|) by using the fact that z/ — g > 7 + p — g = J- □ 

Proceeding along the same line that led to (|5.10p . we obtain a parallel result of 

Lemma 15.71 The proof is postponed to the appendix. 
Lemma 5.13. If 1 < v < |_7 + g\ , then 



tt/2 



ip(y)e(G v {y)) dy 



< 



3tt 



W7 smy u 



(5.21) 



The estimate for the endpoint case v = [ 7 + £> J is essentially the same with the 
argument that led to (|5.21l) . However, the root y v varies with the magnitude of the 
fractional part of 7 + g. Therefore, a more careful treatment is required to obtain a 
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bound that is independent of the magnitude of { 7 + g } . The proof is also postponed 
to the appendix. 

Lemma 5.14. If v = [7 + g J, then 



tt/2 



cp(y)e(G v (y)) dy 



< 47r(iY 7 )- 1 / 2 . 



(5.22) 



Combining these lemmas, we have the following estimate. 
Theorem 5.15. Ift = 0(l), and n + N/2, N/2- 1, then 



\Dy(n,t)\<C 



«2 



Ki + 4k 2 \ V 7 



-(1 + 7) + e(l + £7 + log( 7 + 2)) 



(5.23) 



where C is independent of N,t and 7. 

Proof. Summing up the above three lemmas, we obtain 



L7+eJ+i / 2 



E 



<f{y)e{Gv{y))dy 



L7+eJ-i 
< 2(iV 7 )- 1/2 + Yl 



tt/2 



cp(y)e(G u (y)) dy 



+ 4tt{Nj)- 1/2 + 4(7V 7 )- 2/3 

L7+eJ-i 
< MN!)- 1 ' 2 y — 

J— — ' sin 1 



sm y v 



+ 57r(A^7)" 1/2 + 4(iV 7 )' 2/3 . 



A direct calculation gives 
L7+eJ-i 

sm y v 



L7+0J-1 



e J 1 1 

E — = E 

^ siny^ ^ ,/l - UT- 



^ e) 2 /7 2 
L7+eJ-i +1 ^ 

" S X v/l-(x-p)V7 2 
L7+eJ 1 

: dx 



; dx 



< 



v/i-^-g)^ 2 

7+e 1 



y/l-(x- gf/^ 



dx 



Combining the above two inequalities, we obtain 



L7+eJ+l ^tt/2 



E 

i/=0 



<p(y)e(G v (y))dy 



< 5?r(iV7)- 1/2 (l + 7) + 4(7V 7 )- 2/3 . (5.24) 



It remains to estimate the integral J^ 2 ip(y)e(G v (y)) dy with (p(y) = s'my and 
G v (y) = (iV/7r)(7siny— gy — vy). We choose A = 1/2. The remainder is still bounded 
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by O (e(l + £7 + log(7 + 2))). Obviously, there holds -g - A > -1 since \g\ < 1/2. 
It remains to estimate the shorter sum with v = 0, • ■ ■ ,[j — g + 1/2 \. We deal with 
the cases when v = [j — g \ + 1, v = |_ 7 ~ £ J an d z^ = 0,-- - ,[7 — ,oJ — 1 exactly the 
same with those of Lemmas 15.121 15.141 and 15 . 131 respectively. This results in 



L7-0J+1 r ir/2 



E 

v=0 



(p(y)e(G v (y))dy 



< 57r(7V7)- 1/2 (l + 7) + 4(7V 7 )- 2/3 , 



which together with (|5.24j) gives the final estimate (|5.23j) □ 

6. Estimate of the solution over short time. In this section, we estimate 
the solution over a shorter time interval, i.e., t = 0{e). This is motivated by the 
previous observation that the error already develops to finite magnitude within this 
short time scale. 

Lemma 6.1. If t = 0{e) and n ^ N/2,N/2- 1, then 



\Dy(n,t)\<C 



Kl + 4k 2 



|n+l/2-7V/2|- 2 / 3 + £ 2 / 3 



(6.1) 



where C is independent of N andt. 

The proof of this lemma is essentially the same with that of Theorem 15.151 
Proof. As to f(x) = j/(ne) sm(nx/N) + gx, we have a = g and /? = 7 + g. We 

choose A = 1/2, and the remainder term is bounded by Oie + log(7 + 2)e) = O(e). 

There are only two terms in the shorter sum (|5.3[) . i.e., v = 0, 1 since t — 0(e). When 

v = 0, using Lemma T5. 41 with A = (n + 1/2— L)/ir, we obtain 



Tl/2 



<p{y)e(G (y)) dy 



< 



n+l/2-L 



Using (|5.20[) with v = 1, we obtain 

fTr/2 



tp{y)e(Gx{y))dy 



< 4[7V - (n + 1/2 - L)]- 2/3 < 87V" 2/3 . 



As to f(x) = -f/(ire)sm(irx/N) — gx, we still take A = 1/2, and the remainder 
is still bounded by 0(e). There is only one term in the shorter sum, i.e., v = 0. 
Proceeding along the same line that leads to (|5.20j) . we get, for any 8 £ (0,7r/2), 



r/2 



tp(y)e{G (y)) dy 



< 



2 (n + 1/2 -L)S' 
which is minimized by taking S = 7r 1 / 3 (n + 1/2 — L) -1 / 3 G (0,7r/2). This gives 

y{y)e{G (y))dy < 4(n + 1/2 - Ly 2 '\ 



tt/2 



Summing up the above estimates, we get (|6.1[) . □ 
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We use Euler-MacLaurin formula instead of the truncated Poisson summation 
formula to bound Dy(N/2,t), because this approach gives a more explicit bound 
for the remainder. The starting point is the following first-derivative form of Eulcr- 
MacLaurin formula. For any real valued function f(x) in [a, b] with continuous first 
derivative, we have 



mdx= b —^(f(a) + f(b)) 



2N 



b — a 



N 



I — i \ / 



N 



(x - a)N 
b — a 



fc=l 
(x - a)N 
b — a 



f'(x)dx. 



(6.2) 



Starting with (|5.1|) , and applying Euler-MacLaurin formula (|6.2I) to 
f(x) = sin 2 x cos[(27/e) sinx] 
with a = and b = tt, we obtain 



- g sm - cosKt) = ^ /(*) dx + - ^ — - 



Nx 



The remainder can be directly bounded as 



f 

N 



Nx 

TT 



Nx 



<-(2 + 2 7 /e) 



The integral -f*f (x) dx can be calculated as follows. 



/7T /*' 7r /2 r 7r /2 

/(x) dx = / /(x) dx = / ! 
o Jo Jo 



sm 2 x 



oo 



oo 



o 

2m 



£(-i) r 

m— 
tt/2 



(2m)! 7 



sm 



2m+2 



x dx 



(2m)! 



sin 2 " 1 x dx 



m=0 



1 (2m +1)!! 
(2m)! (2m + 2)!! 

7^ 2m 1 2m + 1 
(m!) 2 2m + 2' 



This implies the following estimate for Dy(N/2,t) when £ = 0(e). 
Lemma 6.2. Ift = 0{e) and n = N/2 or n = N/2 - 1, then 



Dy(n,t) + 



H 2 



2k 2 



ki + 4k 2 hi + 4k 2 
< 27r(e + 7). 



m=0 



7\ 2m 1 2m +1 
el (m!) 2 2m + 2 



(6.3) 



The above estimate means that Dy(N/2,t) is in an 0(e)— neighborhood of a 
constant value that depends on the ratio t/e. 

Remark 6.3. We cannot directly use the above approach to estimate Dy(n,t) 
because an undesirable term 0(n/N) may appear in the bound. 
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7. Discussion. Based on a simple onc-dimcnsional lattice model, we studied 
the effect of ghost forces for dynamic problems. Ghost forces may arise for dynamic 
coupling methods that were derived from energy approximation and Hamilton's prin- 
ciple [U [29l [27l [36] . In our study, based on an approximate model, we show that the 
error develops rather quickly. On the 0(1) time scale, the error is observed in the 
entire domain, and at the interface, the gradient of the error is 0(1). Therefore, the 
influence of the error is more significant than that of the static case. It also suggests 
that the issue of the ghost force may even be more severe than the artificial reflections 
at the interface. The analysis for the original dynamic QC model requires a different 
method, and it will be pursued in our future works. 

Appendix A. Proof of Lemma 15.131 The proof of this lemma is essentially 
the same with that of Lemma 15.71 

Proof For 1 < v < 7, there exists y v G (0,7r/2) such that F v {y v ) = with cosy,, = 
(y — g)/ r y. For any r\ <G (0, min(y„, 7r/2 — y„)) that will be chosen later, we write 



We deal with the three integrals separately Using Lemma I5T41 with A = \F v (y v — rf)\, 
we obtain 





and 



\My,-v)\ 



> 



> 




where we have used Jordan's inequality (|5.7p in the last step. This gives 




Using Lemma 15.41 again with A = \F v {y v +7j)|, we have 
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Furthermore, 



\F v (y v + rf)\ = -F v (y v + n) 



Nj cos y v ~ cos(y, y + 77) 



> 



7r sin(y„ + rj) 
2Nj sin(y„ + 77/2) sin(r?/2) 

7r sin(y„ + 77) 

2^777 sm(y v + 77/2) 
7r 2 sin(y„ + 77) 



which yields 



72 



¥>(z/)e(Gv(j/)) dy 



< 



ir sm(y v + 77) 
2iV777sin(j/ 1/ + ?7/2)' 



If y v G (0,7r/4], wc would require that 77 e (0, y v ). We will have, sin(y„ + 7/) < 
sa\2y v < 2siny t/ , and siny„ < sin(y„ + ?//2) since y„ < y„ + rj/2 < 2y„ < it/2. The 
bound for the above integral is simplified to 



tt/2 



^{y)e{G„{y)) dy 



< 



A trivial bound for the remaining integral yields 

(p{y)e(G v {y)) dy 

-v 

Summing up all the above estimates, we obtain 



./V777 



< 277. 



tt/2 



y{y)e{G v {y)) dy 



< h 277. 

~ 2^777 ' 



(A.l) 



Taking 77 = (Nj) 1 / 2 siny„, which is less than y v , we obtain 



tt/2 



ip{y)e{G v {y)) dy 



< 



3tt 



2 sin y v 



< 



3tt 



2^Nj sin y„ \Nl v ^7 sin y„ 

Now, if yj, € (7r/4,7r/2], we would require that 77 € (0,7r/2 — y v ). We will have 
sin(y„ + 77) > siny„ since y„ < y u +rj/2 < y u +rj < tt/2, and sin(y„+77) < 1 < ^siny^ 
since sin y v > 1/ V2- We bound the second integral as 



tt/2 



<p{y)e(G v {y)) dy 



< 



V2ir 



This yields 



tt/2 



ip(y)e (G„(y)) dy 



< 



.ZV777' 

(1 + V2)tt 
2Njrj 



277. 



In this case, wc might choose 7/ = — (Nj) 1 / 2 provided that 

77 < tt/2 - y v . 
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(A.2) 



With such choice of 77, wc have 



it/2 



(p(y)e(G u (y)) dy 



< 



3tt 

7^7 ' 



This choice of 77 is feasible since 



v - g>l/2> (l/2)y/T/N = 7*7, 

which yields 77 < cosy,, < n/2 — y„, this gives (|A.2[) and completes the proof. 

Appendix B. Proof of Lemma 15. 141 

Proof. There exists y v G (0, 7r/2) such that F v (y v ) = 0, and 

f - 7 - {7 + p} 



cos 



Using the elementary inequality, 



7 



7 



2x x 
1 -<cosx<l xe [0,7r/2]. 

7T 7T 



we obtain 



7r {7 + g} 
2 7 



'tt{7 + 



Proceeding along the same line that leads to (|A.1|) . we get for any 77 G (0, y v ), 



ir/2 



( fi{y)e(Gu(y))dy 



3tt 

: + 2rj. 



We take r\ = (Nj) 1 / 2 , which yields 

»tt/2 



<p(y)e(G u (y))dy 



3tt + 2 4tt 
< — < 



//V7 ~ V]V7 



If { 7 + g } satisfies 



il±*i > 2 (iV7) -i/ 2 , 

7 7T 



then using the left hand side of (|B.2|) . we obtain r\ £ (0, y„). 
On the other hand, if 

7 7T 

then letting S = 2(Nj)~ 1 / 4 , we have 

{ 7 + P } <^ 2 1 — cos s 
~ - 2tt - 2 ' 



7 
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where we have used the right hand side inequality of (|B.1[) . It follows from the above 
inequality and the right hand side of ()B.2I) that y v < 8/ V2 < 8. Using Lemma 15.41 
again with A ~ |-Fi,(<5)|, we get 



tt/2 



ip{y)e{G v {y))Ay 



< 



AF V {8)\ 



We estimate the contribution of the complementary portion of the integral trivially: 



<p{y)e(G u (y)) dy 



< / ydy = —. 



On adding the above estimates we deduce that 



r/2 



<p(y)e(G v {y))dy 



< 



2 n\F v (6)\' 



A direct calculation gives 

Nj cos y v — cos 8 Nj 1 — cos 8 — { 7 + g }/7 



\FAS)\ 
Using (|B.3[) . we obtain 



7r sin 8 



i'm 8 



, „ , „ , A^7 1 - cos (5 iV7 5 7V 7 J 

-PI (S)\ > — ^- = — L tan- > — '—. 

1 Wl - 2tt sin<5 2tt 2 ~ 4tt 



This gives 

r-Tr/2 



<p(y)e{G v (y)) dy 



Finally, we have 

r ir/2 



<p{y)e(G v {y))dy 



< max(47r(7V7)^ 1/2 ,4(7V 7 )- 1/2 ) = 4ir(Nj)- 1/2 
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